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Abstract 

Rayleigh-Benard convection is numerically simulated in two- and three-dimensions 
using a recently developed two-component lattice Boltzmann equation (LBE) 
method. The density field of the second component, which evolves according to the 
advection-diffusion equation of a passive-scalar, is used to simulate the temperature 
field. A body force proportional to the temperature is applied, and the system satis- 
fies the Boussinesq equation except for a slight compressibility. A no-slip, isother- 
mal boundary condition is imposed in the vertical direction, and periodic boundary 
conditions are used in horizontal directions. The critical Rayleigh number for the 
onset of the Rayleigh-Benard convection agrees with the theoretical prediction. As 
the Rayleigh number is increased higher, the steady two-dimensional convection 
rolls become unstable. The wavy instability and aperiodic motion observed, as well 
as the Nusselt number as a function of the Rayleigh number, are in good agreement 
with experimental observations and theoretical predictions. The LBE model is found 
to be efficient, accurate, and numerically stable for the simulation of fluid flows with 
heat and mass transfer. 
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I. INTRODUCTION 



Recently the lattice Boltzmann equation (LBE) method has been developed as a new compu- 
tational fluid dynamics (CFD) method. This method originated form a boolean fluid model known 
as the Lattice Gas Automata (LGA) [|T],^ which simulates the motion of fluids by particles moving 
and colliding on a regular lattice. The averaged fluid variables, such as the density and velocity, 
were shown to satisfy equations similar to the Navier-Stokes equations. The LBE method improves 
this idea by following only the ensemble-averaged distribution functions, therefore eliminating the 
time-consuming statistical average step in the original LGA [^. Simplified collision models were 
later used in place of the collision operator derived from the LGA to improve both the compu- 
tational efficiency and the accuracy. Most noteworthy, the simple collision model of Bhatnagar, 
Gross, and Krook [Q] was applied to the lattice Boltzmann equation, yielding the so called lattice 
BGK model [|],^. The additional flexibility in this approach allows the removal of the artifacts 
of the LGA, specifically the lack of Galilean invariance and the velocity dependent pressure. This 
method was found numerically to be at least as stable, accurate and computationally efficient as 
traditional CFD methods for simulation of simple single-phase incompressible flows More 
importantly, since fluid motion is simulated at the level of the distribution functions, the micro- 
scopic physics of the fluid particles can be incorporate easily like in other particle methods. Many 
complex fluid phenomena due to interparticle interactions, such as capillary phenomena, multiple 
phase flows, and non-linear diffusion, can be simulated naturally [[T0[-|r2|]. 

In most LBE models so far, only mass and momentum conservation is implemented. The 
macroscopic equations of these models correspond to the Navier-Stokes equation with an ideal- 
gas equation of state and a constant temperature. However it is important and sometimes critical 
to have the capability of simulating thermal effects simultaneously with the fluid flows. Obviously 
the temperature distribution in a flow field is of central interest in heat transfer problems. In most 
geophysical flows, the temperature difference is the driving mechanism of the motion of the fluid. 
More importantly, when part of the fluid system undergoes phase transition, as in the boiling and 
evaporation processes, the evolution of the temperature field is directly coupled with the fluid 
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dynamics. Since the LBE method has the most advantage in the simulation of complex fluids 
with multiple phases and phase transitions, it is necessary to develop the capability of simulating 
thermodynamics with the LBE method. 

In general, the simulation of thermal systems by the LBE method has not achieved the same 
success as that of iso-thermal flows. Theoretically, a LBE model with energy conservation can 



be constructed [|13||14|] to yield a temperature evolution equation at the macroscopic level. How- 
ever, the model so obtained suffers from severe numerical instability [[151], especially in three- 
dimensions (3D). Additional stabilization procedures have to be invoked to achieve stability com- 
parable to that of conventional CFD methods, e.g., finite difference schemes. Moreover, when 
interparticle forces are included as in the multi-phase models, the energy conservation is further 
complicated by the potential part of the internal energy. For these reasons constructing a prac- 
tically usable non-ideal-gas LBE model with energy conservation is difficult if not impossible. 
Nevertheless, in many circumstances where the viscous and compressive heating effects can be 
neglected (small Brinkman number limit), the temperature field is passively advected by the fluid 
flow and obeys a much simpler passive-scalar equation. This same equation also governs the diffu- 
sion of each individual component in a fluid mixture. By taking advantage of this formal analogy 
between heat and mass transfer, we can simulate the temperature field as an additional compo- 



nent of the fluid system. Early two-component LGA model [ |16| ] exhibited qualitative features of 
thermal convections. In a previously developed multiple component LBE model [[T7[], we have 
shown that the evolution of the concentration fields is Galilean invariant and obeys Fick's law. 
The diffusivity is independent of the viscosity, allowing a changeable Schmidt number (or Prandtl 
number in the terminology of heat transfer). This model does not implement energy conservation 
and therefore has the same stability as the non-thermal LBE models and other conventional CFD 
methods. By adding one more component, the computation efficiency, either memory-wise or 
time-wise, is not compromised compared with the approach of direct implementation of energy 
conservation because fewer speeds are required for each component. 

In this paper, we present the simulation of the Rayleigh-Benard convection (RBC) as an ex- 
ample. Due to its simplicity and the richness of the phenomena, this problem has been extensively 



studied both theoretically and experimentally [[r^-|22|] and serves as an excellent benchmark prob- 
lem for numerical schemes because detailed results are available for comparison with numerical 
computations. In section |I|, we briefly review the multiple component LBE model and then for- 
mulated it for the simulation of the Boussinesq equation. The implementation of the isothermal 
no-slip boundary condition is also discussed. In section |IJ simulation results are presented and 
compared with theoretical and experimental results. The limitation and some further extensions of 
this method are discussed in section ffVl. 



II. NUMERICAL METHOD 

The following single component lattice Boltzmann equation with BGK collision term describes 
the evolution of the distribution function ^^(x, t) in space x and time t. 

na(x + e„,t + l) -ra„(x,t) = -- [na(x,t) -n^'^'i)(x,t)l , a = l,---,b (1) 

The set of b vectors {ea, a = 1, ■ ■ • , 6} pointing from each lattice site to its neighboring sites 
forms the discretized velocity space of the distribution function. The macroscopic number density, 
n(x, t), and velocity, u(x, t), of the fluid are obtained from Ua as n = Y.a and nu = Y.a ^a^a- 
Eq. (P represents the relaxation of the distribution function to its equilibrium value, n^a^\ which 
is a function of n and u only. The choice of n^^°''> has to ensure that the macroscopic fluid equa- 
tion obtained from Eq. (P by Chapman-Enskog calculation agrees with the Navier-Stokes 



equations. The functional form of n'^^"'^ depends on the structure of the lattice and is usually not 
uniquely determined. For square and cubic lattices in 2D and 3D, the following form of n^^'^ was 
shown to yield Navier-Stokes equations by Qian et al. [Bp: 



ni'"'^ = w„n 



1 + 3e„ ■ u + -[Ba ■ u 



9 , ,o 3u • u 



2' ' 2 



(2) 



Here Wa is a function of le^l and depends on the number of speeds included in the model. In the 
present work, 9 and 15 velocities are used in 2D and 3D computation respectively. The WaS were 
given as 4/9, 1/9, and 1/36 for |e„| = 0, 1, ^2 in 2D and 2/9, 1/9, and 1/72 for [e^l = 0, 1, ^3 



4 



in 3D . It can be easily verify that the 2D distribution function is a degenerate case of the 3D 
version if the flow is two-dimensional. 

A. Multiple component LBE model 

The multiple component LBE model with interparticle interaction [|T^ was originally devel- 
oped for simulation of multi-phase flows and phase transitions. The components can be miscible 
or partially immiscible depending on the strength of the interaction. When the interaction is weak, 
or in a single phase region of a multiphase system, this model can be used to simulate diffusion 
due to various driving mechanisms [|24l]. In this model, the distribution function of each compo- 
nent evolves according to Eq. ([1]). The same form of the equilibrium distribution function given 
by Eq. (0) is used for all the components except that n and u are calculated separately for each 
component. In the absence of any interaction and external forces, the distribution functions of 
all the components were assumed to have a common velocity, u'. The conservation of the total 
momentum at each collision requires that 

s s 

u = 1^ — / L ——^ (3) 

a=l 'o" (T=l 'cr 

where S is the number of components in the system; nia-, t„ and Ua- = J2a the molecular 

mass, the relaxation time, the number density of the component a respectively, and m^n„\i„ = 
m„ Y,a ^a^a is the momentum of component a calculated from its distribution function n^. When 
the force Fa- is applied to component a, the momentum has to be incremented correspondingly. 
This was done by replacing u in Eq. @ with u' + TaF„/ p„. The force F^- in general includes both 
interparticle forces and external forces. For nearest-neighbor interaction, the following form of 
the interparticle force was proposed as it conserves the total momentum of the system and yields 
an adjustable equation of state at the macroscopic level: 

F,^ = -ipa X! 'Sera XI V'<?(x + ea)ea, (4) 

a a 

where is an arbitrary function of the number density of the ath component. 
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In the most general multiple component LBE model with interparticle interaction and external 



nary diffusion, pressure diffusion and forced diffusion. With the equilibrium distribution functions 
given by Eq. (Q), the pressure diffusion does not appear; if a common acceleration is applied to 
all the components, namely Fg. = p^^g, forced diffusion is also absent. The only type of diffusion 
left is the ordinary diffusion due to concentration gradients which obeys Pick's law. In addition, a 
components, e.g. component S, can be made to behave as a passive-scalar by setting its molecular 
mass to zero together with its interaction with all the other components, namely nis and 
Qf^s for a = 1, ■ ■ ■ , S — 1. This component will not contribute to the total momentum of the 
mixture. It is simply advected "passively" and diffuses into the main flow, having no effect on the 
flow. 

For the study of the RBC, we employ a two-component system; component 1 represents the 
motion of the fluid and component 2 simulates a passive temperature field. The distribution func- 
tions of the two components evolve according to Eqs. ([1|) and (g), with u in Eq. being replaced 
by ui + Tcrg for both components. Since the molecular masses of the two components no longer 
appear in the dynamic equations, they are set to unity. The density and the fluid velocity are cal- 
culated from the distribution function of component 1 as p = Z^a'^i ^^id u = ui + g/2, (c.f. 
Ref. JItI]). They satisfy the following equations 



as in the ordinary LBE models. The number density of the second component satisfy the following 
diffusion equation : 



forces, there are three types of diffusions due to different driving mechanisms [^. They are ordi- 




of RBC, it is sufficient to set Qn = 0. The kinematic viscosity u is given by 




(7) 



dn2 



+ V - (nsu) = V - (Wna). 



(8) 



dt 
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The temperature field 9 can be simulated by the density field n2. When the compressibility is 
negligible as in the small Mach number limit, the velocity field is approximately divergence-free 
and the temperature field satisfies the following "passive-scalar" equation: 

BO 

— + u-ve = v-{vve), (9) 

where the diffusivity, V, is given by 

(10) 
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7-2(1 + 9G22i'2dip2/dn2) - ^ 



The diffusivity can be tuned independently of the viscosity by changing either T2 or the interaction 
strength, ^22- For simplicity, Q22 is also set to zero in the present simulation. The LBE model is a 
much simplified version of that in [[1^] since no interparticle interaction is used. 

B. Simulation of the Rayleigh-Benard convection 

la the most common form of RBC, a thin layer of viscous fluid is confined between two 
horizontal rigid boundaries maintained at different temperatures. When the fluid has a positive 
thermal expansion coefficient, and the gravity is in the same direction of the temperature gradient, 
the net buoyancy force is in the opposite direction of the gravity. As the temperature difference 
between the two boundaries is raised above a certain threshold, the static conductive state becomes 
unstable, and convection occurs abruptly. 

The well-known Boussinesq approximation is often used in the study of natural convection. 
With this approximation, the material properties are assumed to be independent of temperature 
except in the body force term, where the fluid density p is assumed to be a linear function of the 
temperature, namely p/poo = I + P{T — Too). Here poo and Too are the density and temperature 
at the reference point, and P the constant thermal expansion coefficient. The gravitational force 
is therefore poog + PooSP{T — Too). After absorbing the first term into the pressure, the effective 
body force is linearly proportional to the temperature variation. 

The Boussinesq equation can be simulated with the two-component LBE model by making 
the external gravitational acceleration, g, in Eq. (H) a linear function of the temperature 9, i.e.. 



g = —qObz, where is the unit vector in the vertical direction and g a parameter controUing the 
strength of the force. In Cartesian coordinates (x, y, z), the lattices are of the sizes x and 
LxXLyX Lz in 2D and 3D respectively. Periodic boundary condition is used in x and y directions, 
and the following no-slip, isothermal boundary condition is used in z direction: 

u = 0,^ = ^ = 0; 

(11) 

u = 0,^ = 1 z^L^. 

Since the LBE fluid is always compressible, an externally applied force will cause a density 
variation. This compressible effect can be eliminated by absorbing into the pressure term the part 
of the body force that corresponding to the body force in the static conductive state, yielding the 
following net external acceleration 

In the conductive state the above external force vanishes and the density field is homogeneous. 

For a given the system size, the characteristic velocity, the Grashof number, the Rayleigh 
number and the Prandtl number are determined by the three parameters ri, T2 and g in the LBE 
model as the following: 

G,^£f, R^GrPr^'§, Pr^^^^. (13) 

The Prandtl number is determined by the two relaxation times used for the two components. Given 
the two basic characteristic dimensionless numbers Pr and R, there is an extra degree of freedom 
in determining t\, ti and g. However, to ensure that the Mach number is small, Vc has to be 
kept small. Once is chosen, all the parameters in the LBE model are determined by the two 
dimensionless numbers R and Pr. 

C. Implementation of the bomidary conditions 

To implement the isothermal, no-slip boundary condition, we must ensure that at the bound- 
ary, the component simulating the fluid flow has zero velocity, and the component simulating 
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the temperature field has fixed density. The mass flux of the second component represents the 
heat transport through the boundaries. Usually the LGA and LBE methods implement the no-slip 
boundary condition by reversing the direction of the incoming particles at the boundary, yielding 
zero averaged velocity. This simple "bounce-back" method was found to be inaccurate [ P5| , P^ . In 
the present work, it results in errors of up to 50% in the critical Rayleigh number. More accurate 
and general methods have been developed to implement velocity boundary conditions in complex 



geometry [[27h29|] . These methods usually involve additional computation at the boundary sites. 
Here, because the boundaries are flat planes, both the velocity and the density boundary conditions 
can be implemented more efficiently. 

When analyzing various implementations of boundary conditions, exact solutions in some sim- 
ple cases are found to be very useful [ p8| , ^ . For simplicity, we consider the time-independent 
one-dimensional situation. All variables depend only on the coordinate perpendicular to the 
wall, so that the spatial dependence can be noted by a single superscript, j, starting from at the 
lower boundary. The elements of the distribution functions can be classified into three groups, 



n\, n 



L and rig, according to the sign of ■ e^. Eq. (|T]) reduces to the following simple form: 



n± — n j_ = — 
r 



7 j(Gq) 
n± — n± 



and nl = ni^'''^\ (14) 



We assume the distribution functions at all sites including the boundary sites are updated uni- 
formly using Eqs. ([T4|). At the lower boundary sites, the groups and Uq are unspecified. The 
only available information about the bulk of the fluid is , from which, n\ is to be constructed 
according to some updating scheme so that certain hydrodynamic boundary condition is satisfied 
at macroscopic level. The "bounce-back" scheme simply sets n]^ = for any a and h satisfying 

= — e^. Obviously this in general does not satisfy Eqs. ( |T4[ ) with u = at the boundary. How- 
ever, if we use the "bounce-back" scheme to calculate the group n^, namely we set n° = for 
any a and h satisfying = — e^, and calculating n\_ using Eqs. with u = and = 6 
in the computation of 7i^}:''°'\ the no-slip boundary condition will be satisfied. Here the summation 
is over all the elements in the group. 

The isothermal boundary condition is imposed by fixing the density of the second component 



at specified values on the boundaries. In the time-independent one-dimensional situation, the 
density profile of the passively convected component can be exactly solved from Eqs. (p^. We 
sum all the elements of the distribution in each of the groups and note the sum as = '^i- By 
summing Eqs. ( [T^ ) we find 

Nt^ -Ni = --\Ni- iVi^'^)" 



T 



(15) 



Using Eq. @ and notice that the velocity only has component parallel to the wall, we can find 
easily that A^^.'-'^'^-' = /Q, independent of the local velocity. Here is the density at the j-th 
position. From the second part of Eq. ((T^), we have + A^i = /3. For any three consecutive j 
values, there are total of seven equations relating the six variables Nj^ to the three density values. 
On eliminating from the seven equations, we find 2n^ = n^~^ + n^^^; namely, the density 
profile is linear in z as the diffusion equation predicts. 
At the boundary sites, we must have 

Nl = N'^/3-N^. (16) 

In the computation, the isothermal boundary condition is implemented by computing the distribu- 
tion function elements in the group n° according to the following equation: 

n° = 2wan - nl (17) 

and then updating it using Eq. ([1]). Here a and b are any pair of indices such that and e;, are 
mirror images of each other with respect to the boundary. In addition to satisfying Eq. ([T^), this 
scheme is also compatible with the no-slip boundary condition. 



m. SIMULATION RESULTS 

We present the simulation results in this section. The two basic characteristic dimensionless 
numbers are the Rayleigh number R and the Prandtl number Pr. The velocity and time reported 
hereafter are in the units of Vc and the characteristic time, Lz/vc, respectively. The diffusion time 
td = Ll/V is VPrR in these units. 
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A. Onset of Rayleigh-Benard instability 



The critical Rayleigh number at which the static conductive state becomes unstable was given 
by the linear stability theory and confirmed by laboratory observations. The static conductive state 
is found to first become unstable to the disturbance of the wave number kc = 3.1 17 in the x-y plane 
when the Rayleigh number exceeds the critical value of = 1707.762. If the deviation from 
the Boussinesq approximation is small, the convection occurs in the form of two-dimensional 
rolls. Since the development of the instability is very slow at near-critical Rayleigh numbers, 
the computation has to be carried out for a long time before stable convection is fully developed. 
Because the first unstable disturbance is two-dimensional, we conduct the near-critical simulations 
primarily in 2D to save CPU time. The results were compared with 3D simulation results for some 
typical cases. 

With periodic boundary condition, the wave number in x-y plane can only take discretized 
values given by 



In 2D, the aspect ratio L^/L^ has to be a multiplication of 2n/kc to accommodate the disturbance 
of the wave number kc. Of course this can only be satisfied approximately on a uniform lattice. 
In the near-critical computation, unless otherwise specified, we chose ~ 2'nLz/kc to save 
computation cost. 

To measure the critical Rayleigh number, computations were started from the static conductive 
state at several different Rayleigh numbers close to Rc- An initial small perturbation was applied 
to the density field. The growth rates of the disturbance were then measured and extrapolated to 
obtain the Rayleigh number corresponding to zero growth rate. 

Shown in Fig. [l| are the typical time-histories of the maximum velocities in z direction for 
three slightly-above-critical Rayleigh numbers of 1720, 1735 and 1750 respectively. The other 
parameters are Pr = 1, ri = 1, and = 50 for all three runs. The peak velocity is found to grow 
exponentially and then saturate at a finite amplitude. The steady-state isotherms and the velocity 




(18) 



11 



field are shown in Fig. ||for the simulation with R = 1750. The growth rates were measured with 
least-square fitting in the exponential growth stage. The fitting results are shown as the straight 
lines. 

The measured growth rates were plotted against the Rayleigh number in Fig. |^. Three sets of 
simulations with ti = 0.55, 1 ., and 1 .5 were performed to investigate the accuracy of using different 
values of ti. All other parameters were the same in these simulations. A solid straight line is fitted 
through the data points for each set of data. The intersections of the lines with the x-axis give the 
Rayleigh numbers corresponding to neutral stability. It is to be seen that near the critical Rayleigh 
number, using a ri other than unity tends to change the growth rates, which causes an error in the 
prediction of the critical Rayleigh number. 

We have also measured the critical Rayleigh number for different Prandtl numbers and with 
different system sizes. The measurement results and the parameters used are summarized in Ta- 
ble |. The biggest error seems to have been caused by using a large r value in the computation. 
Fortunately, this does not impose a significant limitation on the range of physical parameters that 
can be simulated, because for a given Rayleigh number, ti and T2 can always be kept small by 
using a small g. 

Shown in the first part of Table | are five otherwise identical runs with different system sizes. 
The time history of the peak vertical velocity in these runs are plotted in Fig. ^ In the plot, the 
starting times were adjusted so that the initial perturbation levels are the same for all five runs. 
It can be seen that the convergence is fast and differences between the runs with > 20 is 
insignificant. 

Also shown in Fig. ^ is the result of a 3D simulation on a 128 x 128 x 32 lattice with the 
same parameters. The growth rate in the early stage is the same as that in the 2D simulations. 
However, the peak velocity over-shoots before it saturates at the same level. Figs. ^ display a 
series of snapshots of the temperature distribution on a x-y plane laying in the middle between 
the two walls at the times t = 1047, 1320, 1524 and 2273. The grey scales from the darkest to 
the brightest represents the temperature in the range 0.403 < 9 < 0.597. The instability starts in 
the form of an array of convection cells, the superposition of the most unstable mode (k = kc) 
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oriented in the x and y directions, and reaches its maximum near t = 1320 (Fig. The fully 
developed convection rolls oriented in one direction seem to suppress the orthogonal rolls and the 
final convection pattern is purely two-dimensional (Fig. ^). 

B, Higher Rayleigh number 

The two-dimensional convection pattern characterized by the rolls is unstable at higher 
Rayleigh number. As the Rayleigh number is being increased, a series of transitions to more 
complicated states occur, and the form of the convection becomes both three dimensional and 
time-dependent, and eventually turbulent at very high Rayleigh number. Detailed numerical sim- 
ulation of all the complicated transitions and the different forms of convection requires a large 
amount of computation. This is because the form of the convection depends on both the initial 
condition and the boundary conditions. A large number of runs have to be performed to cover the 
parameter space. In addition, the computation has to be carried out for a long time due to the large 
differences among the time scales in the problem. Here we only present the simulation results for 
a few typical situations at moderate Rayleigh numbers due to the limitation of computer resources. 

A two-dimensional simulation at high Rayleigh numbers was performed on a 101 x 50 lattice 
with a Prandtl number of 0.71. The simulation was started from the static conductive state, be- 
ginning with R = 2,000. After the steady-state was reached, the Rayleigh number was raised step 
by step to higher values. The Nusselt numbers measured at the steady states are plotted in Fig. |^ 
against the Rayleigh number. The simulation results of Clever and Busse [|T]] are also plotted 
for comparison. Agreement is found at Rayleigh numbers less than 20,000. At higher Rayleigh 
numbers, the LBE simulation has a lower heat transport. We have raised the Rayleigh number to 
values as high as 10^ for the same resolution. Unlike the thermal LBE model [[T5|], the present 
model remains numerically stable. 

Shown in Figs. ^ are the steady-state isotherms for some typical Rayleigh numbers. As the 
Rayleigh number is increased, the temperature gradient near the boundary becomes sharper; the 
ascending and descending fluid sheets become narrower, and the area in the interior of the fluid 
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with a reversed temperature gradient becomes wider. 

The steady-state solutions were obtained by raising the Rayleigh number gradually after the 
convection roll has established at lower Rayleigh number. It was found however that if the sim- 
ulation is started from the static conductive state with R = 50,000, the system will evolve into 
an oscillatory state. The dominant wavelength is half of that in the steady states solutions shown 
in Figs. ^, and the ascending and descending fluid sheets swing back and forth with a period of 
0.174:Lz/vc. The isotherms at the beginning, the quarter, the half and three quarters of one os- 
cillation period are shown in Fig. ||. This oscillation does not occur in simulations with R < 
30,000. 

Three dimensional simulations were performed for the same Prandtl number and Rayleigh 
numbers on a 128 x 128 x 32 lattice. Again, the computation was started from the static state with 
R = 6,000. Shown in Fig. ^ is the time history of the Nusselt number as the Rayleigh number was 
raised step by step to the values shown on the top of the graph. Greyscale plots of temperature 
distributions on the mid-plane at some typical times for different Rayleigh numbers are shown 
in Figs. The 2D convection rolls have already exhibited some wavy instability at R = 6,000. 
However, the amplitude of the oscillation is so small that the deformation of the convection rolls 
is difficult to be detected from the static plots. To reveal the details of the oscillation, the scale 
has been enlarged and the time history of the Nusselt number replotted in this section. The slow 
decay of the amplitude of the oscillation might be an indication that in agreement with other work- 
ers pT|j32| ], the Rayleigh number of 6,000 is very close to the threshold at which the convection 
becomes time-dependent. The evolution of the convection pattern became more and more compli- 
cated and oscillations of more frequencies were involves as the Rayleigh number was increased. 
At R = 30,000 and 50,000, although the time history of the Nusselt number appeared to be quite 
chaotic, the temperature field plotted in Figs. |1^ still posses rather regular structures and patterns. 
Simulations at higher Rayleigh number would require higher resolution. A detailed investigation 
of the transitions in RBC as the Rayleigh number is increased requires large number of runs and 
is certainly beyond the scope of the present work. 
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IV. CONCLUSIONS AND DISCUSSIONS 



In this paper, we presented a method of simulating temperature evolution in fluid systems 
using multiple component LBE model. By simulating the temperature field using an additional 
component, we are able to avoid the numerical instability plaguing the thermal LBE models. The 
algorithm is simple, and the requirement on computational resources is twice of that for a non- 
thermal LBE code. As an example, we studied the Rayleigh-Benard convection using this method. 
The results agree very well with theoretical predictions and experimental observations both at near- 
critical and moderate Rayleigh numbers. 

The density of the additional component satisfies a passive-scalar equation. In the simulation 
of the Boussinesq equations, the external force is made to be a linear function of this passive 
scalar. However, this passive-scalar can represent other properties of the fluid satisfying more 
complicated equations. More importantly, when the equation of state is coupled with this passive 
scalar, the dynamic process of phase transition can be simulated. We defer the discussion of the 
details to a future publication. 
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TABLES 

TABLE I. Critical Rayleigh number obtained by extrapolating growth rate data at slightly supercritical 
Rayleigh numbers. 
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FIGURES 

FIG. 1. Typical time-histories of the peak vertical velocity in 2D simulation during the onset of the 
instability. The Rayleigh numbers are slightly above i?c- Other parameters are: Pr = 1, ri = 1, and the 
system size is 101 x 50. The sohd straight hnes are the drawn by least-square fitting, the slop of which 
gives the growth rates of the instability. 

FIG. 2. Steady state isotherms and velocity field in a two-dimensional simulation with the Rayleigh 
number of 1750. The resolution is 101 x 50. 

FIG. 3. Growth rates of the instabiUty are found to depend Unearly on the Rayleigh number near Rc. 
The symbols are the results of measurement from the time history of the peak vertical velocity, and the 
straight Unes are fitted through the data points. The exact critical Rayleigh number is obtained by extrapo- 
lating the data to the zero growth rate. Three sets of simulations with different values of ri were performed 
to determine the effects of different ri on the accuracy. The simulations were performed on a 2D 101 x 50 
lattice and is for Pr = 1. 

FIG. 4. Time histories of the peak vertical velocity for different system sizes. It is to be seen that the 
simulation results converge for system sizes > 20. The peak velocity in the 3D simulation has the same 
growth during the early development of the instabiUty and saturates at the same level. The 3D effects peak 
around t = 1320. 

FIG. 5. Grey scale plot of temperature distribution on the mid-plane between the two walls at, (a) t = 
1047, (b) t = 1320, (c) t = 1524 and (d) t = 2273. The gray scales from the darkest to the brightest are 
mapped to the temperature range 0.403 < 9 < 0.597. 

FIG. 6. The steady-state Nusselt number as function of the Rayleigh number in two-D simulations. The 
LBE results agree with that of Clever and Busse for Rayleigh number less than 20,000. 

FIG. 7. Two-D simulation. Isotherms at steady states as the Rayleigh number is raised gradually to (a) 
10,000, (b) 20,000, (c) 30,000 and (d) 50,000. 
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FIG. 8. Isotherms in two-D simulation. The simulation was started from the static conductive state with 
R = 50,000. The system evolves into a oscillatory state. The isotherms are taken at (a) the beginning, (b) 
one quarter, (c) half, and (d) three quarter of one oscillation period. 

FIG. 9. Time history of the Nusselt number in a 3D simulation as the Rayleigh number is increased 
step by step at times indicated by the vertical dashed hnes. For the first segment, the scale has been adjusted 
to the range of 2.08 to 2.082 and the Nusselt number replotted. 

FIG. 10. Greyscale plots of typical temperature distribution on the mid-plane between the two walls for 
Rayleigh number (a) 6,000, (b) 8,000, (c) 10,000, (d) 20,000, (e) 30,000, and (f) 50,000. 
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